Constraining the carbon cycle in JULES-ES-1.0

Code to run analysis and generate figures for McNeall et al. (2023?) “Constraining the carbon cycle in JULES-ES-1.0”

source("JULES-ES-1p0-common-packages.R")
source("JULES-ES-1p0-common-functions.R")
source("JULES-ES-1p0-common-data.R")

0.1 Failure analysis

low_npp_ix <- which(Y[,'npp_nlim_lnd_sum'] < 1e5)
# code from https://stackoverflow.com/questions/28182872/how-to-use-different-sets-of-data-in-lower-and-upper-panel-of-pairs-function-in


#X <- matrix(runif(300), ncol=3)
#Y <- matrix(c(sort(runif(100, 0, 10)), 
#              sort(runif(100, 0, 10)), 
#              sort(runif(100, 0, 10))), ncol=3)

#pdf(file = 'figs/run-failure-pairs.pdf', width = 12, height = 10)
x1 <- X[low_npp_ix, ]
x2 <- X_nlevel0

XY <- rbind(x1, x2)


pairs(XY,
      lower.panel=function(x, y, ...) {
        Xx <- x[seq_len(nrow(x1))] # corresponds to X subset
        Xy <- y[seq_len(nrow(x1))] # corresponds to X subset
        #usr <- par("usr"); on.exit(par(usr))
        #par(usr = c(range(x1[, -ncol(x1)]), range(x1[, -1]))) # set up limits
        points(Xx, Xy, col = zblue, pch = 19, cex = 0.8)
       # if(par('mfg')[2] == 1) axis(2) # if left plot, add left axis
        #if(par('mfg')[1] == ncol(x1)) axis(1) # if bottom plot add bottom axis
      }, 
      upper.panel=function(x, y, ...) {
        Yx <- x[(nrow(x1) + 1):length(x)] # Y subset
        Yy <- y[(nrow(x1) + 1):length(y)] # Y subset
        
        #cntr <- outer(Yx, Yx, FUN='*') # arbitrary function for contour
       # usr <- par("usr"); on.exit(par(usr))
        #par(usr = c(range(x2[, -1]), range(x2[, -ncol(x2)]))) # set up limits
        points(Yx, Yy, col = zred, pch = 19, cex = 0.8)
        #contour(Yx, Yy, cntr, add=TRUE)
        #if(par('mfg')[2] == ncol(x2)) axis(4) # if right plot, add right axis
        #if(par('mfg')[1] == 1) axis(3) # if top plot, add top axis
      }, 
      #tick=FALSE, # suppress the default tick marks
      #line=1,
      gap = 0,
      xlim = c(0,1), ylim = c(0,1),
      labels = 1:d,
      oma = c(2, 18, 2, 2)) # move the default tick labels off the plot 

reset()

legend('left', legend = paste(1:d, colnames(lhs)), cex = 1.1, bty = 'n')
legend('topleft', pch = 19, col = c( zred, zblue), legend = c('failed', 'zero carbon cycle'), bty = 'n', inset = 0.02, cex = 1.1 )

#dev.off()

1 Visualising the ensemble range

It’s important to remember that the design of the experiment is multiplication factors of the original parameters. This might be important for the “hold” value in a sensitivity analysis, as the “standard” value and the median value of the ensemble will not be the same.

#pdf(file = 'figs/lhs_range.pdf', width = 6, height = 8)
par(las = 1, mar = c(5,8,2,1))
lhs_min <- apply(lhs_wave0_wave01_all, 2, min)
lhs_max <- apply(lhs_wave0_wave01_all,2, max)

plot(lhs_max, 1:d, type = 'n', xlim = c(0,10), axes = FALSE, xlab = 'multiplying factor', ylab = '')

abline(v = 0, lty = 'dashed', col = 'grey')

abline(v = 1, lty = 'dashed', col = 'tomato2')

segments(x0 = lhs_min, y0 = 1:d, x1 = lhs_max, y1 = 1:d )

points(lhs_min, 1:d, pch = 20)
points(lhs_max, 1:d, pch = 20)

axis(2, at = 1:d, labels = colnames(lhs))
axis(1)

#dev.off()

1.1 Wave00/Wave01 Ensemble behaviour in key (constraining) outputs.

Global mean for the 20 years at the end of the 20th Century. There is still a significant low bias on cVeg output.

wave00col <- 'skyblue2'
wave01col <- 'tomato2'

wave00col <- 'dodgerblue2'
wave01col <- 'firebrick'
rangecol <- 'grey'

# Histogram of level 1 constraints
hcol = 'darkgrey'
lcol = 'black'

#pdf(file = 'figs/level_2_constraints_hists.pdf', width = 8, height = 8)
par(mfrow = c(2,2), fg = 'darkgrey', las = 1, oma = c(0.1, 0.1, 4, 0.1))

trunc <- function(x, vec){
  
  dat <- x[x < max(vec) & x > min(vec)  ]
  
  dat
  
}


h <- hist(Y_const_level1a_scaled[,'nbp_lnd_sum'], main = 'NBP', xlab = 'GtC/year', col = makeTransparent(wave00col,150))
hist(trunc(Y_const_wave01_scaled [,'nbp_lnd_sum'], h$breaks) ,
     col = makeTransparent(wave01col,150) , breaks = h$breaks, add = TRUE)

rug(Y_const_stan_scaled['nbp_lnd_sum'], lwd = 2)

polygon(x = c(0, 100, 100, 0), y = c(0, 0, 1000, 1000),
        col = makeTransparent(rangecol, 60),
        border = makeTransparent(rangecol))

h <- hist(Y_const_level1a_scaled[,'npp_nlim_lnd_sum'],col = makeTransparent(wave00col,150), main = 'NPP', xlab = 'GtC/year')
hist(trunc(Y_const_wave01_scaled [,'npp_nlim_lnd_sum'], h$breaks) , 
     col = makeTransparent(wave01col) , breaks = h$breaks, add = TRUE)

rug(Y_const_stan_scaled['npp_nlim_lnd_sum'], lwd = 2)

polygon(x = c(35, 80, 80, 35), y = c(0, 0, 1000, 1000),
        col = makeTransparent(rangecol, 60),
        border = makeTransparent(rangecol))

h <- hist(Y_const_level1a_scaled[,'cSoil_lnd_sum'], col = makeTransparent(wave00col,150), main = 'Soil Carbon', xlab = 'GtC')
hist(trunc(Y_const_wave01_scaled [,'cSoil_lnd_sum'], h$breaks) , 
     col = makeTransparent(wave01col,150) , breaks = h$breaks, add = TRUE)

rug(Y_const_stan_scaled['cSoil_lnd_sum'], lwd = 2)

polygon(x = c(750, 3000, 3000, 750), y = c(0, 0, 1000, 1000),
        col = makeTransparent(rangecol, 60),
        border = makeTransparent(rangecol))

h <- hist(Y_const_level1a_scaled[,'cVeg_lnd_sum'], col = makeTransparent(wave00col,150), main = 'Vegetation Carbon', xlab = 'GtC')
hist(trunc(Y_const_wave01_scaled [,'cVeg_lnd_sum'], h$breaks) , 
   col = makeTransparent(wave01col,150)  , breaks = h$breaks, add = TRUE)

rug(Y_const_stan_scaled['cVeg_lnd_sum'], lwd = 2)

polygon(x = c(300, 800, 800, 300), y = c(0, 0, 1000, 1000),
        col = makeTransparent(rangecol, 60),
       border =  makeTransparent(rangecol))

reset()

legend('top', horiz = TRUE, fill = c(makeTransparent(wave00col, 150), makeTransparent(wave01col, 150), makeTransparent(rangecol, 60)), legend = c('Wave00', 'Wave01', 'AW range'))

#dev.off()

1.2 What proportion of wave01 fall within Andy Wiltshire’s constraints?

Just under a third. Points at a significant model discrepency in cVeg

Of the 400 members of the wave01 ensemble, 128 pass Andy Wiltshire’s Level 2 constraints.

length(level2_ix_wave01)

[1] 128

length(level2_ix_wave01) / ntrain_wave01

[1] 0.32

lcol_wave0 <- makeTransparent('dodgerblue2',  120)
lcol_wave01 <- makeTransparent('firebrick',  120)
lcol_wave01_level2 <- 'gold'
stancol = 'black'

linePlotMultiEns <- function(years, ens1, ens2, ens3, col1, col2, col3, ylab, main, ylim = NULL){
  # Plot wave00 and wave01 timeseries on top of one another
  
  nt <- length(years) 
  if(is.null(ylim)){
    
  ylim = range(c(ens1[,1], ens1[,nt], ens2[,1], ens2[ ,nt], ens3[,1], ens3[, nt]))
  }
  
  else ylim <- ylim
  
  matplot(years, t(ens1), type = 'l', lty = 'solid',ylim = ylim, col = col1,
        ylab = ylab, main = main, xlab = '',
        bty = 'n')
  matlines(years, t(ens2), col = col2, lty = 'solid')
    matlines(years, t(ens3), col = col3, lty = 'solid')
}

#pdf(file = 'figs/carbon-cycle-timeseries-waves-constrained.pdf', width = 10, height = 12)
par(mfrow= c(3,5), las = 1, mar = c(4,4,1,0))

linePlotMultiEns(years = years, ens1 = npp_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = npp_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = npp_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'NPP')

lines(years,npp_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 =  nbp_ens_wave00[without_outliers_ix_wave00,], 
                 ens2 = nbp_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = nbp_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01,col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'NBP', ylim = c(-10,10))

lines(years, nbp_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = cSoil_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = cSoil_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = cSoil_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'cSoil', ylim = range(c(cSoil_ens_wave00[,1], cSoil_ens_wave00[,164])))

lines(years, cSoil_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = cVeg_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = cVeg_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = cVeg_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'cVeg')

lines(years, cVeg_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = lai_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = lai_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = lai_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'Lai')

lines(years, lai_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = rh_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = rh_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = rh_lnd_sum_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01,  col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'RH')

lines(years, rh_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = fLuc_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = fLuc_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = fLuc_lnd_sum_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'fLuc')

lines(years, fLuc_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = fHarvest_lnd_sum_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = fHarvest_lnd_sum_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = fHarvest_lnd_sum_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'fHarvest')

lines(years, fHarvest_lnd_sum_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = treeFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = treeFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = treeFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'treefrac'
                 )

lines(years, treeFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = shrubFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = shrubFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = shrubFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'shrubfrac'
)

lines(years, shrubFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = baresoilFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = baresoilFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = baresoilFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'baresoilfrac')

lines(years, baresoilFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, c3PftFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = c3PftFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = c3PftFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'c3PftFrac')

lines(years, c3PftFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, c4PftFrac_lnd_mean_ens_wave00[without_outliers_ix_wave00,],
                 ens2 = c4PftFrac_lnd_mean_ens_wave01[without_outliers_ix_wave01,],
                 ens3 = c4PftFrac_lnd_mean_ens_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'c3PftFrac')

lines(years, c4PftFrac_lnd_mean_stan, col = stancol, lty = 'solid', lwd = 2)

reset()

legend('bottomright', legend = c('wave00','wave01','wave01 level2','standard'), lty = 'solid', lwd = 1.5, col = c(lcol_wave0, lcol_wave01, lcol_wave01_level2, stancol), inset = c(0.05, 0.15) )

#dev.off()

#pdf(file = 'figs/carbon-cycle-timeseries-anomaly-waves-constrained.pdf', width = 10, height = 12)
par(mfrow= c(3,5), las = 1, mar = c(4,4,1,0))

linePlotMultiEns(years = years, ens1 = npp_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = npp_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = npp_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'NPP')

lines(years,npp_stan_anom, col = stancol, lty = 'solid', lwd = 2)
linePlotMultiEns(years = years, ens1 =  nbp_ens_anom_wave00[without_outliers_ix_wave00,], 
                 ens2 = nbp_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = nbp_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01,col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'NBP', ylim = c(-10,10))

lines(years, nbp_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = cSoil_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = cSoil_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = cSoil_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'cSoil', ylim = range(c(cSoil_ens_anom_wave00[,1], cSoil_ens_anom_wave00[,164])))

lines(years, cSoil_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = cVeg_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = cVeg_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = cVeg_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'cVeg')

lines(years, cVeg_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = lai_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = lai_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = lai_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'Lai')

lines(years, lai_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = rh_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = rh_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = rh_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01,  col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'RH')

lines(years, rh_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = fLuc_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = fLuc_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = fLuc_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'fLuc')

lines(years, fLuc_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)
 
linePlotMultiEns(years = years, ens1 = fHarvest_lnd_sum_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = fHarvest_lnd_sum_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = fHarvest_lnd_sum_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = 'GtC', main = 'fHarvest')

lines(years, fHarvest_lnd_sum_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = treeFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = treeFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = treeFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'treefrac'
                 )

lines(years, treeFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = shrubFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = shrubFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = shrubFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'shrubfrac'
)

lines(years, shrubFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)

linePlotMultiEns(years = years, ens1 = baresoilFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = baresoilFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = baresoilFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'baresoilfrac')

lines(years, baresoilFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)


linePlotMultiEns(years = years, c3PftFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = c3PftFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = c3PftFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'c3PftFrac')

lines(years, c3PftFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)


linePlotMultiEns(years = years, c4PftFrac_lnd_mean_ens_anom_wave00[without_outliers_ix_wave00,],
                 ens2 = c4PftFrac_lnd_mean_ens_anom_wave01[without_outliers_ix_wave01,],
                 ens3 = c4PftFrac_lnd_mean_ens_anom_wave01[level2a_ix_wave01, ],
                 col1 = lcol_wave0, col2 = lcol_wave01, col3 = lcol_wave01_level2,
                 ylab = '%', main = 'c3PftFrac')

lines(years, c4PftFrac_lnd_mean_stan_anom, col = stancol, lty = 'solid', lwd = 2)

reset()

legend('bottomright', legend = c('wave00','wave01','wave01 level2','standard'), lty = 'solid', lwd = 1.5, col = c(lcol_wave0, lcol_wave01, lcol_wave01_level2, stancol), inset = c(0.05, 0.15) )


#dev.off()

1.3 Constraining to level 2 with the emulator

nunif <- 50000
X_unif <- samp_unif(nunif, mins = rep(0,32), maxes = rep(1, 32))
colnames(X_unif) <- colnames(X)

# Can this go in common data? Would be needed for checking emulator fits
# Create fit lists for the combined data wave00 level 1a and wave01
Y_const_level1a_wave01_scaled_list <- mat2list(Y_const_level1a_wave01_scaled)

fit_list_const_level1a_wave01 <- mclapply(X = Y_const_level1a_wave01_scaled_list , FUN = km, formula = ~., design = X_level1a_wave01,
                                   mc.cores = 4, control = list(trace = FALSE))

Y_const_level1a_wave01_scaled_pred <- multiPred(Y = Y_const_level1a_wave01_scaled, Xpred = X_unif, fit_list = fit_list_const_level1a_wave01)

level2_ix_em_unif_wave00_wave01 <- which(Y_const_level1a_wave01_scaled_pred$pred_mean[,'nbp_lnd_sum'] > 0 &
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'npp_nlim_lnd_sum'] > 35 & 
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'npp_nlim_lnd_sum'] < 80 &
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'cSoil_lnd_sum'] > 750 &
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'cSoil_lnd_sum'] < 3000 &
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'cVeg_lnd_sum'] > 300 & 
                                           Y_const_level1a_wave01_scaled_pred$pred_mean[,'cVeg_lnd_sum'] < 800
)


(length(level2_ix_em_unif_wave00_wave01) / nunif) * 100

[1] 12.014

X_stan_norm <- normalize(matrix(rep(1, 32), nrow = 1), wrt = lhs)

colnames(X_unif) <- 1:32

#pdf(file = 'figs/pairs_level2_ix_em_unif_wave00_wave01.pdf', width = 12, height = 12)

par(oma = c(0,0,0,3), bg = 'white')

panel_hist_local <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}

pairs(rbind(X_unif[level2_ix_em_unif_wave00_wave01, ], X_stan_norm),
      gap = 0, lower.panel = NULL, xlim = c(0,1), ylim = c(0,1),
      panel = dfunc_up_truth,
      diag.panel = panel_hist_local,
      cex.labels = 1,
      col.axis = 'white',
      dfunc_col = rb
      )


image.plot(legend.only = TRUE,
           zlim = c(0,1),
           col = rb,
           legend.args = list(text = 'Density of model runs matching the criteria', side = 3, line = 1),
           legend.shrink = 0.6,
           horizontal = TRUE
)

legend('left', legend = paste(1:32, colnames(lhs)), cex = 1, bty = 'n')

#dev.off()